'''
Created on Jan 29, 2011

@author: oabalbin
'''

import os
import glob
from optparse import OptionParser

from exome.jobs.base import JOB_SUCCESS, JOB_ERROR
from exome.jobs.job_runner import qsub_cac, qsub_loc, run_local
from exome.jobs.config import ExomePipelineConfig, ExomeAnalysisConfig
from collections import defaultdict, deque

# Global Variables
NODE_MEM=45000.0
NODE_CORES=12
SINGLE_CORE=1
MEM_PER_CORE= int(float(NODE_MEM) / NODE_CORES)
# wt=walltime
WT_SHORT= "24:00:00"
WT_LONG= "60:00:00" #"100:00:00"


def depthOfCoverage(list_input_bam_files, output_bam_file, ref_genome, 
                    interval_list, minimun_coverage, use_mem, path_to_gatk):
    '''
    It produces the read coverage over a 
    particular interval
    
    java -jar /nobackup/med-mctp/sw/bioinfo/gatk/GenomeAnalysisTK-1.0.4905/GenomeAnalysisTK.jar  
    -R /nobackup/med-mctp/sw/alignment_indexes/gatk/hg19/hg19.fa -T DepthOfCoverage 
    -I ../sample/TP53_mut0/TP53_mut0_merged.recal.smdup.bam 
    -I ../sample/TP53_mut1/TP53_mut1_merged.recal.smdup.bam 
    -L chr17:7571720-7590863 -ct 40 -o ./DepthOfCoverage.tmp
    '''

    input_files=[]
    for bamfile in list_input_bam_files:
        input_files += ['-I']+[bamfile]
    
    if interval_list is not None:
        # call in specific sites
        call_on_this_regions=['-L',interval_list]
    else:
        #generate calls in all sites
        call_on_this_regions=['']

        
    

    gatk_command=path_to_gatk+'GenomeAnalysisTK.jar'
    args = ['java','-Xmx'+str(use_mem)+'m', '-jar',gatk_command,'-T',
              'DepthOfCoverage', '-R', ref_genome, '-o', output_bam_file,
              '-ct', minimun_coverage]
    args=args+input_files+call_on_this_regions
    
    args=map(str,args)
    args= [a.replace(',',';') for a in args]
    comd = ",".join(args).replace(',',' ').replace(';',',')

    return comd


def get_depthOfCoverage(analysis, configrun, jobrunfunc, depends, interval_list):
    '''
    '''
    extra_mem, num_cores = configrun.gatk_use_mem, configrun.gatk_num_cores
    path_to_gatk = configrun.gatk_path
    target_exons = configrun.gatk_target_exons
    genomes = configrun.genomes['human']
    ref_genome, snpdb_vcf, indeldb_file = genomes.gatk_ref_genome, genomes.snpdb, \
                                         genomes.indeldb
    my_email=configrun.email_addresses
    
    # to change this for input from the config file
    minimun_coverage=40
                                         
    if  not interval_list:
        interval_list = target_exons
    
    
    cohort_samples = defaultdict(deque)
    for sp in analysis.samples:
        if sp.category=='benign':
            cohort_samples['benign'].append(sp.sorted_mmarkdup_bam)
        else:
            cohort_samples['tumor'].append(sp.sorted_mmarkdup_bam)

    jobn='covg'   
    for sptype, list_bam_samples in cohort_samples.iteritems():
        output_bam_file = os.path.join(analysis.analysis_dir, 'DepthOfCoverage.'+sptype)
        
        comd = depthOfCoverage(list_bam_samples, output_bam_file, ref_genome, 
                        interval_list, minimun_coverage, MEM_PER_CORE, path_to_gatk)
        
        jobidvcf = jobrunfunc(jobn+sptype, comd, SINGLE_CORE, cwd=None, walltime=WT_SHORT, pmem=None, 
                           deps=depends, stdout=None, email_addresses=my_email)
    
    return jobidvcf

        

if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-r", "--config_file", dest="config_file",
                            help="file with run configuration")
    optionparser.add_option("-a", "--analysis_file", dest="analysis_file",
                            help="file with experiment configuration") 
    optionparser.add_option("-l", "--interval_list", dest="interval_list",
                            help="list of intervals to extract the coverage") 

    optionparser.add_option("--local", dest="local", action="store_true", default=False)
    optionparser.add_option("--cluster", dest="cluster", action="store_true", default=False)

    (options, args) = optionparser.parse_args()    

    config = ExomePipelineConfig()
    config.from_xml(options.config_file)
    analysis = ExomeAnalysisConfig()
    analysis.from_xml(options.analysis_file, config.output_dir)
    #Default when called from the command line
    depends=None
    
    if not (options.local ^ options.cluster):
        optionparser.error("Must set either --local or --cluster to run job")
    if options.local:
        jobrunfunc = run_local
    elif options.cluster:
        jobrunfunc = qsub
    
    if not options.interval_list:
        interval_list=[]
    else:
        interval_list=options.interval_list
        
    
    get_depthOfCoverage(analysis, config, jobrunfunc, depends, interval_list)

    
